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We study the g-modes of fast rotating stratified neutron stars in the general relativistic Cowling 
approximation, where we neglect metric perturbations and where the background models take into 
account the buoyant force due to composition gradients. This is the first paper studying this problem 
in a general relativistic framework. In a recent paper I I , a similar study was performed within the 
Newtonian framework, where the authors presented results about the onset of CFS-unstable g-modes 
and the close connection between inertial- and gravity-modes for sufficiently high rotation rates and 
small composition gradients. This correlation arises from the interplay between the buoyant force 
which is the restoring force for g-modes and the Coriolis force which is responsible for the existence 
of inertial modes. In our relativistic treatment of the problem, we find an excellent qualitatively 
agreement with respect to the Newtonian results. 

PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Sf, 97.10. Sj 

I. INTRODUCTION 

Stellar oscillations serve as a unique tool in revealing the nature of the stellar interior and there is a collaborative 
effort during the last three decades to collect this type of information from astrophysically important stars ^2, 
Actually, the Sun's seismic activity has been studied in more detail due to its proximity and thus there already exists an 
independent branch in astrophysics called helioseismology. Since the mid 90s, it has been shown that asteroseismology 
can be combined with gravitational waves in order to shed light in the interior structure of relativistic stars |4l|5l|6j[7]. 
This initial work has been extended by many groups and even the most exotic cases of compact objects like quark 
stars |H] or compact stars in alternative theories of gravity [HI [TD] have been studied. In addition, more detailed 
microphysics has been taken into account, for example, stars with superfluid components [TTJ [12] or even crust and 
magnetic fields [13', [14] . 

Rotation enriches the stellar oscillation spectrum with new frequencies and many complications. This is mainly due 
to the breaking of the spherical symmetry of the problem which complicates the form of the oscillation equations and 
adds to the oscillation spectrum with new families of oscillation modes, such as the so-called inertial modes, while it 
infiuences significantly all other families. In the general relativistic framework, due to the emission of gravitational 
radiation, there appears a new type of instability which is present only in rotating stars, the so-called CFS instability 
[T^ ITS], for a review see [T7]. Thus the rotational modes, r-modes, are unstable to any rotational rate while other 
types of modes like the f- and the g-modes are becoming unstable only for fast rotating neutron stars. Especially, 
the so-called g-modes of rotating stars are present only in stars with internal stratification which may be associated 
with composition variations ^181 and are expected to be important for nascent neutron stars. The g-modes have 
been studied in detail mainly for non-rotating neutron stars in the Newtonian and the general relativistic framework 

[iIlEoJEUES]. 

As it has been mentioned earlier, the effect of rotation on the perturbation equations is dramatic and the majority of 
studies on stellar oscillation and instabilities has been done using either the so-called ' slow rotation approximation ' or 
using Newtonian theory |17[ 123] . The slow rotation approximation assumes that the star is rotating at such a rate that 
the effects of rotation on the shape of the star can be neglected. Only in the last five years or so it became possible to 
study fast rotating stars in the general relativistic framework either using a linearized form of relativistic perturbation 
equations [231 ISS US] or by evolving the non-linear hydrodynamical equations in the relativistic framework [27] [25J [53] . 

The most recent and detailed study of effects of rotation on g-modes by [T] assumes Newtonian gravity but it deals 
properly with the rotation, i.e. the stellar models are ' fast ' rotating. In particular, it has been demonstrated in 
this paper, that the g-modes for high rotation rates become rotationally dominated by the Coriolis force and they 
approach particular inertial modes. In our work, we perform a similar study but in the general relativistic framework 
which describes better both the equilibrium and the perturbed neutron star configurations. Still, one does not expect 
dramatic qualitative changes as one moves from Newton's theory to general relativity, especially for the study of the 
g- and inertial modes and this is what we demonstrate here. Moreover, the Cowling approximation that has been used 
here is known already from the '60s, thanks to the pioneering work of Robe 30J, to be an excellent approximation for 
the study of g-modes and the higher p-modes but it is only qualitatively correct for the study of the lower p-modes 
and the f-mode. 
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In the next two sections we present the equations that describe the problem together with the boundary conditions 
of the problem, while in Section 4 and 5 we review the numerical techniques that have been used and describe the 
results of our simulations. Finally, in the last section we summarize the conclusions and suggest further work. 

II. PROBLEM SETUP 

The basic formulation of this work was already drafted in a previous article by the authors, see |26] (abbreviated 
GK08 in the remainder of this paper) for the case of barotropic oscillations. Here, we will shortly summarize the 
fundamental points and the extensions that were made in order to handle stratified stars. 

In the relativistic Cowling approximation, where one neglects the perturbations of the spacetime, the hydrodynamic 
equations that govern linear oscillations around a static equilibrium take the form 

V^(<5V) = g^''{d^,8T^, - F^^JT^, - ^^JT^x) = , (1) 

where g'"^ is the contravariant metric tensor of the background model, STf^y is the perturbed energy- momentum-tensor 
and F'^K/i are the Christoffel symbols. We further make the assumption that the matter has no viscosity or shear 
stresses. In this case it can be described by a perfect fluid and 6T^i, has the form 

5T^^ = (e +p)(w^(5u^ + u^Su^) + {5p+ 6e)u^u^ + Spg^^, . (2) 

Here e is the energy-density, p is the pressure, is the 4-velocity and Sp, Se, Su^ are its corresponding perturbations. 
Energy density and pressure are not independent from each other but are related via an equation of state (EoS) which 
we assume to be polytropic, i.e. 

p=Kp^+^^^ where e = p+Np. (3) 

Here p is the rest-mass density, K the polytropic constant, N the polytropic exponent and F = 1 + 1/iV the polytropic 
index. 

The system of equations (jlj leads to four relations between the six unknown quantities Su^, 6e and 6p. The 
normalization constraint u^Su^j^ — gives an additional relationship (which takes a very simple form especially in a 
comoving coordinate system, that is 5ut = 0) but still the system is underdetermined. What is needed is an auxiliary 
correlation between the Eulerian pressure perturbation 5p and the energy density perturbation 6e. For adiabatic 
oscillations and a polytropic equation of state for both the background and the perturbation this relation is given by 
(see [31] for example) 

Here Fi is the polytropic index of the perturbed fluid elements and <^^, are the two relevant components of the 
fluid displacement vector in cylindrical coordinates. Since the stationary background is axisymmetric, there is no 
pressure variation with respect to (p and hence no additional term proportional to in equation d4| . In the barotropic 
case when F = Fi, the square of the speed of sound is given by = Fp/(e + p) and relation M reduces to the well 
known form 6p = c^Se also used in GK08. 

We now have a connection between dp and Se but at the cost of introducing two new quantities, i.e. the fluid 
displacement components, that have to be evolved in time as well. There has to be a relation between Jit^ and fluid 
displacement since the latter is basically an integrated velocity. In the Cowling approximation this dependency is 
given by (see e.g. [52] ) 

Su^ = (.9^^ + Uf_,u^)LuC , (5) 

where denotes the Lie derivative along m^. This Lie derivative will introduce a time derivative on the righthand 
side of equation ^ which then will serve as an evolution equation for the components of the fluid displacement vector 

III. EVOLUTION EQUATIONS AND BOUNDARY CONDITIONS 

As in GK08 we are using a comoving frame of reference in cylindrical coordinates {g, (, (p, t). In such a system the 
metric reads 

ds^ = e-^'^ [e^^ (dg^ + dC^) + W^d(p^] - e^^ {dt + adpf , (6) 



3 



where U, k, W and a are axisymmetric metric potentials that depend on g and C only. The (</?, t)-component of the 
metric tensor is proportional to the potential a and vanishes in the absence of rotation. 

With the polytropic equation of state in the form of relation ^ , the polytropic index T is equivalent to the adiabatic 
index and it is 

p ^ e + p dp 
p de 

On the other hand, for stratified systems the speed of sound Cs is defined by 

2 Ti dp 



Using the last two equations one finally arrives at 



r de 



2 

c„ = — — 



e + p 

for stratified stars. This allows us to rewrite equation Q into 

5p = cl5e + clgi + clg2 , (7) 

where we definded 



Since both relations Q, ([t]) are linear and equation ([t]) reduces to the barotropic case for gi = g2 = 0, we expect to 
see additional terms in the resulting set of time-evolution equations which are similar to the corresponding expressions 
for the energy density perturbation 5e already obtained in the barotropic version of the problem. This will also serve 
as a simple check of our calculations. 

For fast rotating configurations, it is not reasonable any more to decompose the fluid perturbations into spherical 
harmonics. Centrifugal forces, though only an C'(n^)-efFect where il is the rotation rate of the star, will lead to a 
flattening of the rapidly spinning object where spherical harmonics are no longer appropriate as a system of base 
functions. Instead, we separate the azimuthal dependence into complex exponential functions and we write 

{e+p)We''5u^ = /l(^>,C,^)e™^ 

(e + p)W^e^fec = /2(^>,C,^)e™^ (9) 
{e + p)5u^ = /3(^^,C,^)e™'^ 

for the standard fluid perturbation variables and 

c^e^5i - gi{gX,t)e'^'' (10) 

for the displacement vector variables. Inserting ([t]) into ([2|, using equations (|9|, ( |10| ) and the normalization constraint 
Sut = for replacing ST^^i, in relation ([TJ, finally leads to the system of evolution equations for the fluid perturbations 
in stratified stars: 



dt \ dg dg dg J W dg cj. dg 



dt \dC dC dC J W dC c2 dC 

dh , 2iu. T..2... ^^ Wacle'^'J-^'' fdh df2\ e^--^\^da ^ e''--^\da ^ ,,,, 

-J2 



ie+p)F ie + p)F 

dH im , 2 w f 2 4C/..r , _ , ^^ , Wcle^^-'^^ fdfi , df2\ cy^-^^ a da ^ c^e^^-^fe „ 



{e + p)FW •'^ {e + p)FW ' 
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where 



F := a^cle^^ ~ 



is just a convenient abbreviation whereas the coefhents 



are closely related to the vector field 



1 1 

f " f\ 



Vlnp (12) 



whose amplitude is referred to as Schwarzschild discriminant and which determines the convective stability of g-modes. 
For |A| = we are again in the barotropic limit while |A| < means convective stability and |A| > signals unstable 
g-modes (see e.g. [33]). 

A couple of observations can be made here. Firstly, if we restrict the problem to barotropic perturbations again, then 



gi=g2 = A = B = Q and the evolution equations (111 reduce to the system of equations successfully implemented 



already in GK08. Secondly, the appearance of terms like H + gi + g2 in the last two equations of system (111 as 
well as its derivatives with respect to g and C in the first two equations is a direct consequence of the linearity of the 
perturbation equations and the decomposition made in equation 



The system (111 is still incomplete though. On the righthand side, the fluid displacement varibles gi and 52 appear 
and one needs evolution equations for them as well in order to update these quantities for the next time step. This 
is done with the help of relation ([s]). Since the Lie derivative of a function agrees with the standard differentiation 
of a scalar variable along a vector field and since u'^ has a very simple form in comoving cylindrical coordinates, i.e. 
u'^ = (0, 0, 0, e"*^) , one immediately arrives at the evolution equations for the fluid displacement quantities which 
take the form 

dg2 _ Se^t^ 
~dt ~ {e + p)w'^' 

As already discussed in Section ITT) these equations basically describe that the fluid displacement is obtained by 



integrating the corresponding velocities. Actually, relationship ( 13 1 has already been used for deriving the evolution 
equations for /s and H, since on their righthand side additional time derivatives of gi and 32 appeared; i.e. the 
coefficients proportional to A, B. 

The system of evolution equations is closed by the appropriate boundary conditions that prescribe the behaviour 
of the perturbations at the boundaries of the numerical domain. Using the same computational setup as described in 
GK08, these boundaries consist of the rotation axis and the surface of the star. The original boundary conditions for 
the barotropic evolution variables do not change at all in the case of stratification and a very similar analysis can be 
carried out for the additional quantities gi, 52 • Let us first consider the surface of the star. Restricting ourselves to 
linear perturbations, both variables gi and 52 in equation (|7| have to be finite in order to add up to a finite Eulerian 



pressure perturbation 5p at the surface. By virtue of equation ( 10 ) this means, that gi and 32 have to vanish there. 
The boundary condition for all evolution quantities along the surface therefore is 

yi|surfacc y2|surfacc /slsuifacc -^|surfacc (-^^) 

^?l|surfacc 52|surfacc ■ 

For the rotation axis, one has to distinguish between scalar and vectorial perturbations the same way it was done 
in the barotropic case. Scalar perturbations have to be unique along the axis for all values of m when varying the 
azimuthal angle Lp while vectorial perturbations are only allowed to change like cos{ip) or sin((y5) near the axis. Put in 
other words, only for m = ±1 they can attain nonzero values there. In the case presented here, gi which represents 
clearly has vectorial character while (72 can be treated as a scalar quantity concerning the proper behaviour along 
the rotation axis. Table |l] gives a complete overview of the boundary conditions for all evolution quantities. 



IV. NUMERICAL IMPLEMENTATION 



The implementation to solve the system of equations ( 11 ) is virtually adopted from GK08. Just like in that paper, 
a pseudospectral code is used to generate the axisymmetric background model which then gets interpolated onto a 
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TABLE I: Boundary conditions for the perturbation quantities along the rotation axis. Here, 'f&c' stands for 'finite and 
continuous'. 



m- value /i|axis 


/'2I axis 


/"si axis 


axis 9l 1 axis 


<72|axis 











f&c 


f&c 


±1 





f&c 


f&c 





else 















TABLE II: Background model parameters; is the rotation rate, M the gravitational mass and Vp/re the ratio of polar to 
equatorial circumferential radius 

model Q, (kHz) Vp/ve M (Mq) re (km) 



BUO 


0.0 


1.0 


1.4 


14, 


,16 


BUOl 


0.694 


0.995 


1.4 


14, 


,19 


BU02 


0.981 


0.99 


1.41 


14, 


,22 


BU03 


1.695 


0.97 


1.42 


14, 


,36 


BUI 


2.182 


0.95 


1.43 


14, 


,52 


BU2 


3.062 


0.9 


1.47 


14, 


,92 


BUS 


3.712 


0.85 


1.5 


15, 


,39 


BU4 


4.229 


0.8 


1.54 


15, 


,91 


BUS 


4.647 


0.75 


1.59 


16, 


,53 


BU6 


4.976 


0.7 


1.63 


17, 


,27 


BU7 


5.214 


0.65 


1.67 


18, 


,16 


BUS 


5.344 


0.6 


1.69 


19, 


,29 



finite difference grid that extends the original computational domain in polar direction to < < tt. The time 
evolution equations themselves are solved by an Iterated Crank-Nicholson scheme with swapped weights (^H [55] ) 
and a weighting factor of a = 0.6. Similar to the barotropic case, the code is numerically unstable. Exponentially 
growing modes destroy the stability already after very few timesteps. A dissipative, second order Kreiss-Oliger term 
was added to the fluid displacement variables gi and 32 as well. However, it was found that the amount of dissipation 
needed for gi, g2 is at least two orders of magnitude smaller than for the rest of the evolution variables. With this 
artificial viscosity, the numerical evolution stays stable for all simulations shown in this paper. The typical evolution 
time was around 150 ms; this leads to a frequency resolution of A/ « 7 Hz. 

In addition, mode recycling routines as in [57], [5S] and GK08 were used to identify and enhance specific modes as 
well as for the extraction of eigenfunctions. 

V. RESULTS 
A. Background Models and Initial Data 

In this paper, we use the polytropic BU model series ([35]) for which K = 100, F = 2 and a fixed central rest-mass 
density of pc = 1-28 x 10~^ in units of G = c = Mq = 1. In the convential cgs-system this means a central rest-mass 
density of pc — 0.79 x 10^^ g/cm^ and a corresponding energy density of Ec — 0.89 x 10^^ g/cm^. In the nonrotating 
limit this leads to a neutron star with 'standard values' for gravitational mass, M = IAAIq, and circumferential 
radius, r — 14.16 km, and is labelled BUO. Starting from there, the angular velocity is subsequently increased in 
such a way that the ratio of polar to equatorial radius rp/r^ decreases by a factor of 0.05 for every new equilibrium 
configuration until the mass-shedding limit of f2 = 5.36 kHz is reached. Table [H] gives an overview over the model 
parameters for this specific equation of state. 

In order to resolve the splitting of the co- and counterrotating g-mode oscillations accurately at low frequencies, 
additional equilibrium configurations were taken into account between the nonrotating model and BUI. Since the 
latter one already has a spin frequency of around 40% of the Kepler- limit, three auxiliary equilibrium models were 
constructed to trace the regime of low rotation rates. 

Concerning initial data, we are following [T] and mostly use a Gaussian radial profile in the pressure perturbation 
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of the form 

<5p = Apexpfl^) YU9,(t>). (15) 

Ws{d)J 

Here, A is a dimensionless amplitude factor, p the rest-mass density, rs{9) the coordinate radius of the star at polar 
angle 9 and rg, q determine the centre and width of the Gaussian perturbation. This type of initial data typically 
excites not only g-modes but also pressure- and inertial modes. Applying mode-recycling techniques is crucial to 
extract, identify and follow specific modes throughout the range of possible spin frequencies. 



B. Barotropic Oscillations 

A thorough test of the numerical implementation of system ( [TT| in the barotropic case, i.e. F = Fi, can already 
be found in GK08 since the time evolution equations for stratified stars derived in this paper reduce to the form 
studied in GK08 for barotropic oscillations (see equation (11) there). Figure [l] shows some typical results of such 
a simulation. The left panel depicts the time evolution of a perturbed velocity component for the most rapidly 
rotating stellar model BUS. Keep in mind, that we have to use artificial viscosity in order to stabilize the code against 
spurious modes. This also allows for very long integration times; in the case of Figure [l] the simulation was cancelled 
after 0.25 seconds although longer runs pose no problems. The artificial viscosity has a secondary effect though. In 
addition to the numerical dissipation that is always present when using finite-difference representations of differential 
equations, the extra viscosity amplifies this dissipative effect which can clearly be seen in the decreasing amplitude 
of the initial perturbation. In frequency space, this typically leads to a broadening of the less dominant peaks in 
the power spectrum. However, subsequent mode recycling runs for the specific mode of interest can countervail this 
behaviour. 
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FIG. 1: Left panel: Time evolution of the perturbed p-component of the 4- velocity using the fastest spinning background model 
BUS running for about 250 milliseconds. Right panel: Mode splitting for two fundamental, nonaxisymmetric pressure modes as 
well as the change in frequency for the purely axial Z = m = 2 inertial mode; all frequencies are shown in the comoving frame 
of reference. 



C. Stratified Oscillations 



Once the perturbations of the fluid obey a different equation of state than the unperturbed background, buoyancy 
can drive a new class of oscillation modes. It acts in a similar manner as usual convection, i.e. perturbed fluid elements 
rise in the star as long as their densities are smaller than the unperturbed surroundings. As soon as the densities 
are equal, there is no net force acting on the fluid elements but due to their inertia, they overshoot this equilbrium 
position and move into regions where they are heavier than the surroundings. They begin to fall, overshoot again but 
now into regions where they are lighter and the cycle starts again. 

It is already well known from Newtonian calculations (see e.g. |33j). that p- and g-modes differ significantly both 
in their frequencies and their eigenfunctions. While the former have frequencies at around 1.5 kHz and higher for 
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TABLE III: Comparison of the first / = 2, 3, 4 fundamental g-modes with the results from Yoshida & Lee (2000). The frequencies 
are given in units of (j/\/Gpc, where pc is the central rest-mass density. The literature values are denoted with YL. 



1 


51 91 




YL 


2 


0.0189 0.0181 


3 


0.0227 0.0223 


4 


0.0255 0.0254 



typical equations of state and are virtually confined in a region between an inner reflection point and the surface of 
the star, g-modes have much lower frequencies, depending on the degree of stratification, and the large amplitudes 
that they reach in the stellar interior decrease pretty strong towards the neutron stars' surface. Figure [2] shows a 
comparison between the two-dimensional eigenfunctions of the fundamental p- and g-mode, reconstructed from our 
relativistic time-evolution code and a composition gradient of Fi = 2.05. 




FIG. 2: Comparison between pressure- and gravity-mode eigenfunctions. Left panel: Energy-density eigenfunction for the 
I = A,m = 2 fundamental p-mode. Right panel: Energy-density eigenfunction for the Z = 4, m = 2 fundamental g-mode. 



In the case depicted in Figure |2] the simulation was performed on the non- rotating BUO model with trial initial data 
that excited the ''^f-mode at a frequency of 2.6 kHz as well as the ^gi-mode at 304 Hz. The layout of the computational 
domain and the physical meaning of its coordinates (ct, r) were already described in GK08. Essentially cr is a radial 
coordinate ranging from the origin of the star to its surface while r is an angular coordinate that reaches from one 
pole of the star, i.e. = Q in. a spherical system, to the other one at 9 — tt, where t — 1 corresponds to the equatorial 
plane. As one can see clearly, the maximum amplitude of the pressure mode is indeed located near the surface of the 
star while the g-mode eigenfunction has its peak value near the center of the star and falls off quickly towards the 
surface. 

In order to further check the reliability of our code with respect to these new non-barotropic oscillations only 
found in stratified stars, we computed the frequencies of the fundamental Z = 2, 3, 4 g-mode for a nonrotating, weakly 
stratified star with F = 2.0 and Fi = 2.0004. A similar calculation has been performed in [S7] within the Newtonian 
framework. The following Table HI summarizes the results. The overall agreement is very good, although it should be 
noted that the small value of Fi leads to correspondingly low frequencies for the fundamental g-modes. Together with 
the constant rest-mass density of the BU model series of pc = 0.79 x 10^^ g/cm^, this means that the absolute values of 
the three g-mode frequencies lies within a bandwith of w 20 — 30 Hz. It is very demanding to accurately determine 
the frequency peaks in this small frequency range. High order modes even have smaller frequencies; according to [37], 
the first overtones for Z = 2, 3, 4 can be found in the range i' ~ 14 — 20 Hz. Currently, the numerical noise in our code 
prevents us from properly determining these low frequencies. 

In [T], the fundamental I — 2 g-mode as well as several overtones were computed for different values of Fi. These 
values go as high as Fi = 2.4 and lead to considerably higher frequencies for stratified oscillations. We calculated 
the mode frequencies for the fundamental and the first overtone of the I = 2 g-mode of the non-rotating and non- 
barotropic BUO model for which F = 2.0 and checked it again with the literature values. The results are depicted in 
Table IV and show very good agreement; only for larger composition gradients, i.e. larger Fi, there are some larger 
deviations. Actually, the authors in [1] used a relativistic, time-independent eigenvalue formulation for slowly rotating 
stars in the Cowling approximation developed and described in |38] and a new Newtonian time-evolution approach 
for rapidly rotating configurations to show the good agreement between these two independent methods in the case 
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TABLE IV: Comparison of the fundamental and first overtone 1 = 2 g-mode with the Passamonti et al. (2009) results. The 
frequencies are given in units of a/^/Gp2, where pc is the central rest- mass density. The literature values from the boundary 
value approach are denoted with PHAJH. Additionally, the absolute values of the mode frequencies (in Hz) for the model BUO 
are shown in parentheses. 



Pi 


'91 


^51 


^52 


'92 




PHAJH 




PHAJH 




2.05 


0.195 


0.196 (227) 


0.134 


0.132 (153) 


2.2 


0.383 


0.382 (442) 


0.266 


0.258 (299) 


2.4 


0.531 


0.503 (582) 


0.370 


0.344 (398) 



of a nonrotating star. The relative difference between these two schemes never exceeds 7%. The time evolutions 
performed with our code agree a little bit better with the eigenvalue results, which should be no surprise since for 
the special case of a nonrotating star, both relativistic codes are solving the same problem; either as boundary value 
problem as in [38] or as a time-dependent evolution problem presented in this paper. 

To summarize, the results compared so far for nonrotating stars and different values of the composition gradient are 
consistent with previous studies within the error bounds. We can now turn to the investigation of rotational effects 
on the frequencies of g-modes and study the influence of stratification on these frequencies. 

From the nonrotating and slow-rotation limit it is already well known that the angular part of the fluid perturbations 
and, in the case of a coupled spacetime evolution, the corresponding metric perturbations, can be decomposed in terms 
of scalar, vector and tensor spherical harmonics which split into two distinct classes depending on their behaviour 
under space reflection, i.e. {r,9,(f) {r,Tr — 9,tt + (p). Polar perturbations that are expressed in terms of Yim and 
its gradient VYzm transform as (—1)' while axial perturbations behave like r x VYim and transform with (— 1)'+^. In 
a non-rotating, barotropic, perfect fluid star there are no non-trivial, i.e. zero-frequency, axial modes. In the case of 
the Cowling approximation, the only modes present are the pressure modes which have polar parity and frequencies 
of around 1 kHz and higher (see Figure [T]). This changes, once rotation or a composition gradient is included. Both 
cases will introduce a new class of modes; inertial modes which are restored by the Coriolis force and the g-modes 
where gravity is the restoring force due to buoyancy. 

In contrast to inertial modes which are degenerate for a nonrotating background star, gravity modes have a non- 
vanishing oscillation frequency already in the nonrotating limit. The following Figure [3] shows a comparison between 
fluid perturbations of a barotropic BUO-model and a stratified BUO star with a composition gradient of Fi = 2.1. 
The two upper panels show an 0.02 sees excerpt from a time-evolution performed for these different configurations. 

In both cases one can count roughly 37 oscillations of the dominant mode during this time which leads to an 
estimation of the corresponding frequency of about 1.85 kHz and in fact can be identified with the fundamental 
quadrupolar mode (see also Figure [T|) . But in the stratified case depicted in the right part of Figure |3] there is an 
additional low-frequency modulation at around 300 Hz which is not present in the barotropic time-evolution and 
which is the fundamental quadrupolar g-mode at 314 Hz. 

Both mode classes, gravity- and inertial modes, have mixed polar and axial components and both of them occupy 
the low- frequency regime (see Figure [l] for the purely axial I — m = 2 inertial mode and Table 111 IV for g-modes) . 
Figure |4] shows the variation of the fundamental g-mode frequency for several values of the composition gradient Fi. 
The left panel shows the polar-led ^gi-mode while the right panel depicts the axial-led ^gi-mode. Both panels show 
the oscillation frequencies in the inertial frame and for comparison, the barotropic fundamental pressure mode is also 
displayed for the corresponding harmonic indices. In both cases, the g-modes behave in a similar fashion. The higher 
the composition gradient Fi, the higher are also the mode frequencies; this can clearly be seen from Figure [4] 

In contrast, the pressure modes only show a negligible effect on the magnitude of stratification, although their 
frequencies are also slightly larger. For example, in the nonrotating case, the barotropic ^f-mode has a frequency of 
V = 1.879 kHz. For Fi = 2.05 this changes to = 1.88 kHz and is absolutely neglectable while for Fi — 2.3 the 
frequency shifts to = 1.886 kHz which is still only a relative shift of 0.4%. This picture gets modified though a 
little bit in the case of rotation, mostly because the counterrotating mode moves down in frequency and therefore 
the relative difference increases. However, this difference never exceeds the 7%-level. The following Figure [5] shows 
a comparison between the frequency variation of the fundamental quadrupolar ^f-mode with rotation for a purely 
barotropic star and a stratified background model with a comparatively large composition gradient of Fi = 2.3. As 
already stated, the effect of stratification on pressure modes are minute. 

As soon as the frequency of a specific mode drops down to zero, i.e. when it becomes degenerate, this oscillation mode 
is prone to the so-called CFS-instability. This type of rotational dragging instability was originally discovered in [13], 
thoroughly investigated in [161 139] and always works when there is a coupling mechanism to some radiation field (for 
example gravitational radiation) and a rotating background that leads to a splitting between co- and counterrotating 
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f/Hz f/Hz 



FIG. 3: Time evolution of fluid perturbations for a barotropic and a stratified equilibrium model in the nonrotating limit. Top 
left panel: Oscillation of 5ug for a barotropic neutron star with F = Fi = 2.0. Top right panel: Oscillation of 5ut^ for a stratified 
neutron star with F = 2.0 and Fi = 2.1; the modulation in the velocity perturbation is due to the influence of the fundamental 
^gi-mode. Bottom left panel: The power spectral density for the barotropic model. Bottom right panel: The corresponding 
PSD plot for the stratified star. 

modes. A rotational instability can only occur when the mode frequency vanishes along a sequence of equilibrium 
models. 

It has already been observed in GK08, that the ^f-mode of the model BU series goes marginally unstable just at the 
mass-shedding limit. This can also be seen from the left panel of Figure |4] and Figure [5] where the dotted horizontal 
line represents the zero- frequency limit. The quadrupolar fundamental mode is a very good gravitational wave emitter, 
however it does not enter the instability unless the background star spins extremely fast. The absolute values for 
the non-axisymmetric ^f-mode will change in a fully relativistic treatment of the problem with coupled spacetime 
evolution (see for example |28p. but qualitatively the picture remains unchanged; the CFS- unstable quadrupolar 
f-mode usually only works in rapidly rotating stars, if at all. 

On the other hand, due to their low frequency, the fundamental g-mode becomes unstable at low rotation rates 
as can be clearly seen in Figure |4j Depending on the degree of stratification, the critical spin frequency increases 
for higher composition gradients since the mode frequencies are also growing in this case. In |40j . Lai gives a spin 
frequency estimation for the onset of the secular g-mode instability depending on the nonrotating oscillation frequency 
(To which is 

n^ = QMao, (16) 

where ctq = 2tivo. For the six unstable modes in Figure |4] we did least square fittings with all available data points 
and a quadratic fitting polynomial. We find a somewhat smaller value than 0.68 in |40] which is f2c = cgq with 
c — 0.58 ± 0.01. This is roughly 17% smaller than the original value and agrees very well with Newtonian results 
which lead to a similar decrease. However, it should also be pointed out that the g-mode instability is not thought 
of being a very significant gravitational wave emitter; viscosity will damp away the instability very soon. Actually, 
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FIG. 4: Variation of mode frequencies for non-barotropic rotating stellar models with F = 2.0 and several values of Fi. Left 
panel: Frequencies for the I = m — 2 barotropic f-mode and fundamental g-mode for different composition gradients. Right 
panel: Frequencies for the I = 3,m — 2 barotropic f-mode and fundamental g-mode for different composition gradients; all 
frequencies are shown in the inertial frame of reference. 




the growth times for unstable g-modes, in the slow rotation approximation and by using a post-Newtonian approach, 
have already been calculated in |2J- The results there suggest growth times of the order of 10^ - 10^ sees, which is 
considerably longer than the time viscosity will need to suppress completely the instability. 

For stratified stars, there is now the interesting situation, that two classes of modes practically occupy the same 
domain in frequency space. Both emerge from degenerate, zero frequency oscillations and as the rotation rate of the 
background star increases, both classes will be affected by the growing Coriolis force. Inertial modes are inherently 
dominated by the Coriolis force but of course this force will also influence the g-modes. For these type of oscillation, 
buoyancy is the restoring force but the faster the compact object is spinning, the more important the Coriolis force 
will become. How dominant it gets essentially depends on two parameters. First of all there is the rotation rate; the 
Coriolis force will grow with increasing angular velocity of the background star. One would therefore expect that in 
the case of rapidly spinning background models, the g-modes will become more similar to pure barotropic inertial 
modes and this in fact verified by Newtonian calculations; see again [1]. However, the rotation rate and hence the 
strength of the Coriolis force is limited by the mass-shedding limit; it cannot grow beyond the breakup frequency. 
This is where the second parameter, which determines the interplay beween g-modes and inertial modes, becomes 
important; it is the degree of stratification. In the case discussed in this paper, where polytropic equations of state are 
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FIG. 6: Variation of mode frequencies with the rotation rate of the background star for the counterrotating ^gi-mode and 
several composition gradients Fi. Additionally, a purely barotropic inertial mode for which Fi = F is also depicted and shows 
how g- and inertial mode frequencies approach for high rotation rates and small composition gradients; all frequencies are 
shown in a comoving frame of reference. 



used for both the static equilibrium configuration and the perturbations, it is the difference between Fi and F that 
is crucial. The larger this difference is the larger is also the density difference between unperturbed and perturbed 
fluid and therefore the larger is the effect of buoyancy. This means that if we just choose the polytropic index of 
the perturbed fluid Fi high enough, then the influence of the increasing Coriolis force will be counterbalanced and 
possibly overtaken by the buoyant force. In this case, the similarities between g-modes and barotropic inertial modes 
concerning both frequencies and eigenfunctions may not be so pronounced as it is the case for smaller composition 
gradients. 

We checked this expected behaviour with the unstable, counterrotating fundamental ^gi-mode and various degrees 
of stratification. Since in this case the mode is travelling in opposite direction with respect to the rotation of the star, 
it easily gets dragged forward at a rather low spin frequency and is therefore affected by the CFS-instability. In the 
left panel of Figure [4] these unstable g-modes are depicted in an inertial frame of reference and we are focusing now 
on the unstable branches at low frequencies. Figure [6] shows the result of the simulations for composition gradients 
ranging from Fi = 2.0 for the barotropic inertial mode to Fi — 2.4 for the gravity mode with the highest degree of 
stratification. In contrast to Figure [4] all frequencies are depicted in a reference frame corotating with the star. Figure 
[6] essentially confirms our assumptions about the interplay between g- and inertial modes for growing rotation rates 
and increasing Coriolis force opposed to the growing influence of the buoyant force for larger composition gradients. 
As already discussed earlier, the g-mode frequencies increase for larger degrees of stratification and they also increase 
for higher rotation rates. For small composition gradients, i.e. with Fi = 2.05 and Fi — 2.1 in Figure [6] the mass- 
shedding limit for our chosen background models is high enough for the Coriolis force to dominate over buoyancy in 
the rapid rotation regime. Consequently, the g-modes become similar to certain barotropic inertial modes which in 
particular means that their oscillation frequencies approach each other. The Coriolis force is also relevant for g-modes 
with a higher degree of stratification but in these cases, i.e. with Fi roughly starting from a value of 2.2 in Figure 
[6] the effect of buoyancy is strong enough even at the Kepler-limit to counteract the Coriolis force. This leads to a 
somwhat larger deviation of the g-mode frequencies from the purely barotropic inertial mode frequencies as can be 
clearly seen in Figure [6j 

The similarity between inertial and g-modes for a sufficiently large Coriolis force does not only affect the oscillation 
frequencies but also the mode eigenfunctions. In fact, this is how we are able to find the barotropic inertial mode that 
closestly approaches the Fi = 2.05 counterrotating ^gi-mode. For this, one typically starts at a rather high rotation 
rate of the background model and tries to excite the fundamental g-mode. Several mode-recycling runs are performed 
to obtain an accurate frequency determination. Then the eigenfunctions for the various perturbation variables are 
extracted and put back as initial data but now for barotropic oscillations with Fi = F. This will usually excite some 
inertial modes within the expected frequency range which then again have to be identified by their eigenfunction. 
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FIG. 7: Eigenfunctions of the (5itg-perturbation for the counterrotating ■^pi-mode and the barotropic inertial mode at several 
rotation rates. First row: The eigenfunction of the gravity mode for model BUO. Second row, left panel: The eigenfunction of 
the fundamental gravity mode for model BUOl. Second row, right panel: The corresponding eigenfunction of the barotropic 
inertial mode for model BUOl. Last row, left panel: The eigenfunction of the fundamental gravity mode for model BU6. Last 
row, right panel: The eigenfunction of the inertial mode for model BU6. 



Once the correct oscillation mode has been picked out, one can follow its frequency variation towards lower and higher 
rotation rates. 

Figure [7] shows the results of such simulations for certain eigenfunctions at several rotation rates. The depiction is 
similar to the eigenfunction plots from Figure [2] but this time the amplitude of the eigenfunction is projected onto 
the two-dimensional computational domain rather than additionally plotted along a third axis. Instead, contour lines 
are used to underline the characteristics of the various eigenfunctions. The coordinates on the numerical domain are 
identical to the ones used in Figure [2] i.e. a represents the radial coordinate and ranges from to 1 while r is the 
the angular coordinate and ranges from to 2. The topmost row shows the (5ug-eigenfunction of the fundamental 
counterrotating ^gi-mode for Fi ~ 2.05 and a nonrotating background model at a frequency of Vg = 227 Hz. Since 
there are no inertial modes for 17 = 0, there is also no corresponding eigenfunction. The two panels in the second row 
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FIG. 8: The frequencies of the co- and counterrotating ^gi-mode as well as the barotropic inertial modes ii, 12 which approach 
the gravity modes for high rotation rates as seen from a comoving coordinate system. Similar to polar-led g-modes, the 
increasing Coriolis force dominates over buoyancy for small composition gradients. In contrast, the purely axial r-mode shown 
here is essentially unaffected by the degree of stratification. 



show again the eigenfunction for 6ug but this time for model BUOl rotating at = 694 Hz. The left picture depicts 
once more the fundamental quadrupolar g-mode now at i^g — 251 Hz where one can already see the formation of two 
smaller secondary peaks located symmetrically above and below the equatorial plane. The right panel in the same row 
shows the corresponding eigenfunction of the barotropic inertial mode at the same rotation rate and with a frequency 
of Ui = 111 Hz. Here, the secondary peaks are much more pronounced and almost equal in amplitude compared to 
the main peak. One can already figure out a certain global resemblance between the inertial and g-mode; at least the 
major characteristics are similar in both cases. This similarity becomes even more evident at very high rotation rates. 
The last two panels in the bottom row depict the very same eigenfunctions but now at a spin frequency of f2 = 4.976 
kHz; this is roughly 93% of the mass-shedding limit. In accordance with the frequencies that approach each other for 
high rotation rates, see Figure [6] the same happens with the eigenfunctions. In the last two panels, the ^gi-mode has 
a frequency of i/g — 827 Hz while the barotropic inertial mode oscillates at Ui — 805 Hz. In both cases, the main peak 
gets narrowed in a confined region along the equatorial plane while the secondary peaks grow in amplitude for the 
g-mode and decrease for the inertial mode. By this gradual transition, the eigenfunctions resemble each other very 
good at high spin frequencies. 

A similar observation can be made for axial-led g-modes like the '^gi-mode, whose frequency variation can be seen in 
the right panel of Figure |4] Figure [8] shows both the counter- and corotating branches of this mode for Fi = 2.05 but 
now in a coordinate system comoving with the star. As already described in GK08, the high frequency branch can be 
identified with the retrograde rotating oscillations, i.e. these modes become unstable once the rotation rates exceeds 
a certain threshold, while the low frequency branch corresponds to modes that travel in the prograde direction. What 
is also depicted in Figure [8] are the two barotropic inertial modes that approach both of these branches as the spin 
frequency increases. As for polar-led g-modes, the rather small degree of stratification leads to a dominance of the 
Coriolis force over buoyancy at sufficiently large rotation rates and therefore the g-modes approach eigenfrequencies 
of certain barotropic inertial modes. In addition, we also show in this picture the sensitivity of the purely axial 
I = 2,171 — 2 r-mode frequencies on several composition gradients. As already observed in [1], there is virtually no 
effect of the degree of stratification on the r-mode. Even for the rather large value of Fi — 2.3 the change in frequency 
is essentially negligible. 



VI. CONCLUSIONS 



In this work it as been presented, for the first time, a study of the g-modes for fast rotating stars in the general 
relativistic framework. The results are in qualitative agreement with those by [1] which were derived also for fast 
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rotating stars but in the Newtonian framework. 

This study demonstrates that potentially we are able to analyze future observational data from neutron stars in the 
general relativistic framework for equilibrium models which where constructed by using Einstein's theory of gravity. 

One of the main results observed in this work is the mixing of the g-modes with inertial modes for high rotation 
rates which suggests that it will be difficult for asteroseismology to discriminate the effects of stratification. On the 
other hand, for lower rotation rates the degree of stratification may become observable since the two families of modes 
maintain their distinct character. 

As a next step we plan to study the effect of differential rotation on all families of fluid modes. It is expected 
that differential rotation will be present during the very first moments of the formation of neutron stars when their 
fluid actually will be stratified. In this case the neutron stars are still quite warm, i.e. temperatures considerably 
higher than lO^K while the matter is in the form of degenerate protons p, neutrons n and electrons e which are not 
yet superfluid. This will be an extension of this work and the work that has been done for differentially and slowly 
rotating stars in [38j. 
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